/*Copyright (C) 2010-2012 RAAFCOM B.V. All rights reserved.

This file is part of Velsto and may be used under the terms of:
The GNU General Public License version 3 or higher, as published by
the Free Software Foundation.

This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.*/
/*********************************************************************************** 
 *
 *                     ** VELSIM **
 *
 * Gaussian simulation of time-depth conversion in layer-cake model
 *
 *
 * v1.0		05/01/2012 M.A.van Ravesteijn(Raafcom B.V.)
 *			Initial version
 * v1.1     14/05/2012 M.A.van Ravesteijn(Raafcom B.V.)
 *          added Lineair Regression function
 * v1.2     16/10/2012 M.A.van Ravesteijn(Raafcom B.V.)
 *          char length projectname and resultfilename 300
 *
 * Format Input file (stdin)
 * -- Project name
 * -- Grids file name (Grid format x,y,twt1,twt2....twtnrlayers, other grid parameters 
 *    defined in following rows)
 * -- Gridvvn Grids file array index (layer 1..10) (value -1 no grid parameter)
 * -- Gridvvndev Grids file array index (layer 1..10) (value -1 no grid parameter)
 * -- Gridtwtdev Grids file array index (layer 1..10) (value -1 no grid parameter)
 * -- Avz(Average depth) results filename
 * -- Sdz(Standard Deviation depth) results filename
 * -- Report filename
 * -- nrlayers
 * -- Simulation option
 *    (1 = Single pont single run
 *     2 = Single point simulation
 *     3 = All points single run.
 *     4 = all points full simulation)
 * -- Nr of realizations in simulations(2-2000)
 * -- Nr of layers velocity model.
 * -- (Per layer)
 *   -- Layernr.
 *   -- Velocity type(V0/Vint) 1 = Vint, 2 = Reg Vint, 3 = V0 and k, 4 = Reg V0 and k, 5 = Lin. Regr.
 *                                                                                         (dZ=F.dT + G).
 *   -- Velocity(V0/Vint)(m/s)/F(m/s).
 *   -- Standard Deviation Velocity(m/s)/F(m/s).
 *   -- K(1/s)/G(m)
 *   -- Standard deviation TWT(ms)
 ***********************************************************************************
 */

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <errno.h>
#include <math.h>
#include <time.h>

#define CONS_V			1
#define REG_V			2
#define CONS_V_K		3
#define REG_V_CONS_K	4
#define LIN_REGR        5

FILE *fd_rep_vs,*fd_grids,*fd_avz_vs,*fd_sdz_vs;
FILE *fp_log_vs;

double random_vs();
double invnor_vs(double);
double calc_deltaz(int,double,double,double,double,double);
double round(double);


int main()
{
   double *z1_vs,*z2_vs,*vn_vs,*v1_vs,*v2_vs,*t1_vs,*t2_vs;

   float  velo_vs[10],sdvn_vs[10],rka_vs[10];
   float  sdtwt_vs[10],twt_vs[10],temp2_vs[40],sdtwt_rep[10];
   float  temp3_vs[22],rav_vs[12],rsd_vs[12],regv_vs[10];
   float  xcoord_vs,ycoord_vs,xcoordv_vs,ycoordv_vs;
   float  fvalue_vs;
	
   double sum_vs,sum2_vs,deltaz_vs,deltat_vs;
   double diff_vs,dumt_vs,xt_vs,qt_vs,dumv_vs,xv_vs,qv_vs,tsq_vs,ssq_vs;

   int    ivm_vs[20];
   int    ln_vs[20];
   int    isimul_vs,nrsimul_vs,np_vs,nlay_vs,nsim_vs;
   int	  nrgrid;
   int    igridv0[10],igridvvndev[10],igridtwtdev[10];
   int    i,j,k,l,i1,v;
   int    iloop_vs; 
   int    nv_vs;
   int    nfirst_vs;

   char   projnam[301],fgrids[301],favz_vs[301],fsdz_vs[301],frep_vs[301];
   char   buf[256];
   char   string[500],string2[20];
   char   string3[20];
   char   *p,*q,*end;
   char   *starttime,*endtime;

   typedef int Boolean;
   enum { 	False = 0,
		    True  = 1};
   Boolean skip_point,gridvvn[10],gridvvndev[10],gridtwtdev[10];

   skip_point = False;   
   nv_vs = 10;        
   for (v=0;v<nv_vs;v++) {
        gridvvn[v]=False;
        gridvvndev[v]=False;
		gridtwtdev[v]=False;
   }
    srand(RAND_MAX/2);
	nlay_vs   = 0;
	isimul_vs = 0;
	xcoord_vs  = 0.0;
	ycoord_vs  = 0.0;
	xcoordv_vs = 0.0;
	ycoordv_vs = 0.0;
    fp_log_vs = fopen("velsim.log","w");
	if(	(z1_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(z2_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(vn_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(v1_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(v2_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(t1_vs = (double *) calloc(2000,sizeof(double))) == NULL ||
		(t2_vs = (double *) calloc(2000,sizeof(double))) == NULL 
	   )
	{
		fprintf(fp_log_vs,"Not possible to allocate memory for arrays\n");
		exit(1);
	}
   
   fprintf(fp_log_vs,"In velsim.c\n");
   char startdateStr [9];
   char starttimeStr [9];
   _strdate( startdateStr);
   fprintf(fp_log_vs, "The current date is %s \n", startdateStr);
   _strtime_s( starttimeStr );
   fprintf(fp_log_vs, "The current time is %s \n", starttimeStr);
//  Read Projectname
   gets(projnam);
   projnam[strlen(projnam)]='\0';	
   fprintf(fp_log_vs,"projnam [%s]\n",projnam);
//  Read filename of Grids file
   gets(fgrids);
   fgrids[strlen(fgrids)]='\0';        
   fprintf(fp_log_vs,"fgrids[%s]\n",fgrids);
   fprintf(fp_log_vs,"nv_vs[%d]\n",nv_vs);
//  Read gridparameters
	for (v=0;v<nv_vs;v++) {
	   scanf("%d",&igridtwtdev[v]);
   }
   for (v=0;v<nv_vs;v++) {
	   scanf("%d",&igridv0[v]);
   }
   for (v=0;v<nv_vs;v++) {
	   scanf("%d",&igridvvndev[v]);
   }
/*
 * Read output filenames
 */
   scanf("\n");
   gets(favz_vs);
   favz_vs[strlen(favz_vs)]='\0';	
   fprintf(fp_log_vs,"favz_vs = %s\n",favz_vs);
     
   scanf("\n");
   gets(fsdz_vs);
   fsdz_vs[strlen(fsdz_vs)]='\0';	
   fprintf(fp_log_vs,"fsdz_vs = %s\n",fsdz_vs);
          
   gets(frep_vs);
   frep_vs[strlen(frep_vs)]='\0';	
   fprintf(fp_log_vs,"frep_vs = %s\n",frep_vs);     
      
   if((fd_rep_vs=fopen(frep_vs,"w")) == NULL) {
      fprintf(fp_log_vs,"**ERROR** REP file %s Open failed\n",frep_vs);
	  exit(1);
   }
//   Read nrlayers,simulationoption,nr.realisations,layer parameters
   scanf("%d",&nlay_vs);
   if (nlay_vs > 10) {
	  nlay_vs=10;
	  fprintf(fp_log_vs," Nr.layers>10 set to 10\n");
   }   
	scanf("%d",&isimul_vs);
	scanf("%d",&nrsimul_vs);
	for (i=0;i<nlay_vs;i++) {
	   scanf("%3d%1d %f %f %f %f",&ln_vs[i],&ivm_vs[i],
		  &velo_vs[i],&sdvn_vs[i],&rka_vs[i],&sdtwt_vs[i]);
//   convert standard deviation twt from ms to s 
	   sdtwt_rep[i] = sdtwt_vs[i];
	   sdtwt_vs[i] = sdtwt_vs[i]/1000.0;
	}
 
	fprintf(fd_rep_vs,"************************************\n");
	fprintf(fd_rep_vs,"*** Report on VELSIM Run ***********\n");
	fprintf(fd_rep_vs,"************************************\n");
	fprintf(fd_rep_vs,"\n");
	fprintf(fd_rep_vs,"*** Input ***\n");
	fprintf(fd_rep_vs,"Project         : %s\n",projnam);
    fprintf(fd_rep_vs,"Start Simulation %s %s\n", startdateStr, starttimeStr);
    fprintf(fd_rep_vs,"Nr of Layers    : %2d\n",nlay_vs);
	switch(isimul_vs)
    {
	     case 1:
			fprintf(fd_rep_vs,"Simulation Type : Single Point Single Run\n");
	        break;
	     case 2:
			fprintf(fd_rep_vs,"Simulation Type : Single Point Simulation\n");
	        break;
	     case 3:
			fprintf(fd_rep_vs,"Simulation Type : All Points Single Run\n");
	        break;
	     case 4:
			fprintf(fd_rep_vs,"Simulation Type : Full Simulation\n");
	        break;
	     default: 
			 fprintf(fp_log_vs,"**ERROR** Unknown simulation type %d\n",isimul_vs);
	        exit(1);
	     break;
	}
	fprintf(fd_rep_vs,"Nunber of Realisations : %4d\n",nrsimul_vs);
	fprintf(fd_rep_vs,"\n");
	fprintf(fd_rep_vs,"Layer Parameters\n");
	fprintf(fd_rep_vs,
    " ---------------------------------------------------------- \n"
	" |LY| V |V0/INT|  F   |STDV V|STDV F|  K   |   G   | STDV   |\n"
	" |NR|TYP|(m/s) |(m/s) |(m/s) |(m/s) |(1/s) |  (m)  | TWT(ms)|\n"
	"  -- --- ------ -----  ------ ------ ------ ------- -------  \n");
	for (i=0;i<nlay_vs;i++) {
		if (ivm_vs[i]== LIN_REGR) {
			fprintf(fd_rep_vs," |%2d|%3d|      |%6.1f|      |%6.1f|      |%6.4f|%6.4f|\n",
			ln_vs[i],ivm_vs[i],velo_vs[i],sdvn_vs[i],rka_vs[i],
			sdtwt_rep[i]);
		}
		else {
            fprintf(fd_rep_vs," |%2d|%3d|%6.1f|      |%6.1f|      |%6.4f|       |%6.4f|\n",
			ln_vs[i],ivm_vs[i],velo_vs[i],sdvn_vs[i],rka_vs[i],
			sdtwt_rep[i]);
		}
	}
	fprintf(fd_rep_vs,"\n");
//   Open Grids file
	if((fd_grids=fopen(fgrids,"rb")) == NULL) {
	   fprintf(fd_rep_vs,"**ERROR** Grids file %s Open failed\n",fgrids);
	   exit(1);
	}
//   Report Velocity types
	for (i=0;i<nlay_vs;i++) {
	   i1=i+1;
	   switch(ivm_vs[i])
           {
	     case CONS_V:
	        fprintf(fd_rep_vs,"\n Layer %2d has a Constant VINT",i1);
	        break;
	     case REG_V:
	        fprintf(fd_rep_vs,"\n Layer %2d has a Regionally varying VINT",i1);
	        break;
	     case CONS_V_K:
	        fprintf(fd_rep_vs,"\n Layer %2d has Constant V0 and K",i1);
	        break;
	     case REG_V_CONS_K:
	        fprintf(fd_rep_vs,"\n Layer %2d has Regionally varying V0 and Constant K",i1);
	        break;
		 case LIN_REGR:
	        fprintf(fd_rep_vs,"\n Layer %2d has Lineair Regression function (dZ=F.dT + G)",i1);
	        break;
	     default: 
			 fprintf(fp_log_vs,"**ERROR** Unknown velocity type %d for Layer %2d : %d\n",ivm_vs[i],i1);
	        exit(1);
	        break;
	   }
   	}

	fprintf(fd_rep_vs,"\n");
//   Check if Grid VVN is needed
	for (i=0;i<nlay_vs;i++) {
	   if( ivm_vs[i]== REG_V || 
	       ivm_vs[i]== REG_V_CONS_K )
	   { 
	     gridvvn[i] = True;
	   }
	}
//   Check if VVN Dev Grid file is needed
    fprintf(fd_rep_vs,"\n");
	for (i=0;i<nlay_vs;i++) {
	   i1=i+1;
	   if( igridvvndev[i] > -1 ) {
	      fprintf(fd_rep_vs,"Layer %2d has Regionally varying V0/VINT Uncertainty\n",i1);
 	      gridvvndev[i] = True;
 	   }
 	   else {
 	      fprintf(fd_rep_vs,"Layer %2d has Constant V0/VINT Uncertainty\n",i1);
 	   }
	}
//   Check if TWT Dev Grid file is needed
    fprintf(fd_rep_vs,"\n");
	for (i=0;i<nlay_vs;i++) {
	   i1=i+1;
	   if(igridtwtdev[i] > -1) {
	      fprintf(fd_rep_vs,"Layer %2d has Regionally varying TWT Uncertainty\n",i1);
 	      gridtwtdev[i] = True;
 	   }
 	   else {
 	      fprintf(fd_rep_vs,"Layer %2d has Constant TWT Uncertainty\n",i1);
 	   }
	}

	if((fd_avz_vs=fopen(favz_vs,"w")) == NULL) {
	   fprintf(fd_rep_vs,"**ERROR** AVZ output file %s Open failed\n",favz_vs);
	   exit(1);
	}

	if((fd_sdz_vs=fopen(fsdz_vs,"w")) == NULL) {
	   fprintf(fd_rep_vs,"**ERROR** SDZ output file %s Open failed\n",fsdz_vs);
	   exit(1);
	}
//   set number of realisations following simulation type

	if (isimul_vs == 1 || isimul_vs == 3) {
	   nsim_vs = 1;
	}
	else if (isimul_vs == 2 || isimul_vs == 4) {
	   nsim_vs = nrsimul_vs;
	}
//   set standard deviations to 0 if single run
	if (nsim_vs == 1) {
	   for (i=0;i<nlay_vs;i++) {
	      sdvn_vs[i] = 0.0;
	      sdtwt_vs[i]= 0.0;
	   }
	}

	fclose(fd_grids);

	if((fd_grids=fopen(fgrids,"rb")) == NULL) {
	   fprintf(fd_rep_vs,"**ERROR** Grids file %s Open failed\n",fgrids);
	   exit(1);
	}

	
	fprintf(fd_rep_vs,"\n");
//   Begin while loop over xy points
	np_vs =0;
	nfirst_vs = 0;
    while(fgets(string,500,fd_grids) !=0) {
	   skip_point=False;
       np_vs++;
	   sscanf(string,"%f %f",&xcoord_vs,&ycoord_vs);
       p    = string;
       end  = string + strlen(string) - 1;
       i    = 0;
       while (*p != '\n')
       {  q    = end;
          fvalue_vs = (float) strtod(p, &q);
          if (strcmp(p,q)==0){
             *p='\n';
             break;
          }
          if (i > 1) {
             temp2_vs[i-2] = fvalue_vs ;                    
          }
          i++;
          p    = q;
	   }
//      Set xy output arrays
	   rav_vs[0] = xcoord_vs;
	   rav_vs[1] = ycoord_vs;
	   rsd_vs[0] = xcoord_vs;
	   rsd_vs[1] = ycoord_vs;

	   for (i=0;i<nlay_vs;i++) {
	      if (temp2_vs[i]<1.0 || temp2_vs[i]>9998.0) {
		     fprintf(fd_rep_vs,"Point X= %10.1f Y= %10.1f"
			 " skipped because of invalid TWT's\n",xcoord_vs,ycoord_vs);
	         skip_point = True;
			 if (isimul_vs==1 || isimul_vs==2) {
				 exit(1);
			 }
	      }
	   }
//      Check integrity TWT's of Point
       for (i=1;i<nlay_vs;i++) {
	      diff_vs = temp2_vs[i] - temp2_vs[i-1];
	      if (diff_vs>-2.0 && diff_vs<0.0) temp2_vs[i]=temp2_vs[i-1];
	      if (diff_vs<-2.0) {
		 fprintf(fd_rep_vs,"Point X= %10.1f Y= %10.1f"
			 " skipped because of Discrepancy in TWT's between Layers\n",xcoord_vs,ycoord_vs);
		 skip_point = True;
			 if (isimul_vs==1 || isimul_vs==2) {
				 exit(1);
			 }
	      }
	   }
//	    Get Regional velocities from Grids file
	   for (i=0;i<nlay_vs;i++) {
	      temp3_vs[i] = 0.0;
       }
       for (v=0;v<nv_vs;v++) {      
	      if (gridvvn[v]) {
			    temp3_vs[v] = temp2_vs[nlay_vs + igridv0[v]];
          }
       }
	   for (i=0;i<nlay_vs;i++) {
		  if (gridvvn[i]) {
			  regv_vs[i] = temp3_vs[i];
			  if (regv_vs[i] < 1.0 || regv_vs[i] > 9998.0) {
                 fprintf(fd_rep_vs,"Point X=%10.1f Y=%10.1f skipped because of invalid Reg.V\n",
					 xcoord_vs,ycoord_vs);
		         skip_point=True;
		      }
	      }
	   }
//      Read Velocity Uncertainty Grid Value from Grids file record    
	   for (i=0;i<nlay_vs;i++) {
	      temp3_vs[i] = 0.0;
       }
       for (v=0;v<nv_vs;v++) {      
	      if (gridvvndev[v]) {
     		 temp3_vs[v] = temp2_vs[nlay_vs + igridvvndev[v]];
          }
       }

	   for (i=0;i<nlay_vs;i++) {
		  if (gridvvndev[i] && nsim_vs>1) {
			 sdvn_vs[i] = temp3_vs[i];
		  }
		  else if (nsim_vs==1) {
			 sdvn_vs[i]=0;
          }
       }
//      Read TWT uncertaintie from Grids file           
	   for (i=0;i<nlay_vs;i++) {
	      temp3_vs[i] = 0.0;
       }
       for (v=0;v<nv_vs;v++) {      
	      if (gridtwtdev[v]) {
			  temp3_vs[v] = temp2_vs[nlay_vs + igridtwtdev[v]];
          }
       }
	   for (i=0;i<nlay_vs;i++) {
		  if (gridtwtdev[i] && nsim_vs>1) {
//         Convert sd twt from ms to s
			 sdtwt_vs[i] = temp3_vs[i]/1000.0;
          }
	      else if (nsim_vs==1) {
			 sdtwt_vs[i]=0.0;
	  	  }
      }
		 

	   if (!skip_point) {
	      for (i=0;i<nlay_vs;i++) {              
//          Convert twt's from ms to s
		         twt_vs[i] = temp2_vs[i]/1000.0;
              
          }


	      for (j=0;j<nsim_vs;j++) {
		     v1_vs[j] = 0.0;
		     v2_vs[j] = 0.0;
		     t1_vs[j] = 0.0;
		     t2_vs[j] = 0.0;
		     z1_vs[j] = 0.0;
		     z2_vs[j] = 0.0;
	      }
//         Start loop over Layers
	  for (i=0;i<nlay_vs;i++) {
//         Copy points from previous layer
         for (j=0;j<nsim_vs;j++) {
	         t1_vs[j] = t2_vs[j];
	         v1_vs[j] = v2_vs[j];
	         z1_vs[j] = z2_vs[j];
	     }

//        Get deltat
		 if (i==0) {
		    deltat_vs = twt_vs[0];
		 }
		 else {
		    deltat_vs = twt_vs[i] - twt_vs[i-1];
		 }
//        If deltat < 2*sdtwt do not simulate
		 if (i>0 && deltat_vs <= (2*sdtwt_vs[i])) {
		    rsd_vs[i+2] = rsd_vs[i+1];
		    if (twt_vs[i] == twt_vs[i-1]) {
		        rav_vs[i+2] = rav_vs[i+1];
		    }
		    else {
			   deltaz_vs = calc_deltaz(ivm_vs[i],deltat_vs,velo_vs[i],regv_vs[i],
					                    rka_vs[i],rav_vs[i+1]);

		       rav_vs[i+2] = rav_vs[i+1] + deltaz_vs;
		       for (j=0;j<nsim_vs;j++) {
			      t2_vs[j] = t1_vs[j] + deltat_vs;
			      v2_vs[j] = (2*deltaz_vs)/deltat_vs;
			      z2_vs[j] = z1_vs[j]+deltaz_vs;
		       }
		    }
		 }
//        Else simulate
		 else {

//           Fill time draws
		    if (sdtwt_vs[i] > 0.0001) {
		       if (i > 0) {
				  for (k=0;k<nsim_vs;k++) {
				     diff_vs = -999.9;
				     iloop_vs = 0;
				     while (diff_vs <= 0.0 && iloop_vs < 500) {
                        iloop_vs++;
					    qt_vs   = random_vs();
					    xt_vs   = invnor_vs(qt_vs);
					    dumt_vs = twt_vs[i] + xt_vs*sdtwt_vs[i];
					    diff_vs = dumt_vs - t1_vs[k];
			         }
			         t2_vs[k] = dumt_vs;
                     if (diff_vs <= 0) t2_vs[k] = t1_vs[k] + 0.0001;
			      }
		       }
		       else {
			      for (k=0;k<nsim_vs;k++) {
			         qt_vs   = random_vs();
			         xt_vs   = invnor_vs(qt_vs);
			         t2_vs[k] = twt_vs[i] + xt_vs*sdtwt_vs[i];
				  }
		       }
		    }
		    else {
		       for (k=0;k<nsim_vs;k++) {
			      t2_vs[k] = twt_vs[i];
		       }
		    }

//           Fill velocity draws
			double vel = 0.0;
			if( ivm_vs[i]== REG_V || 
	            ivm_vs[i]== REG_V_CONS_K )
            {
			   vel = regv_vs[i];
			}
	        else
			{
		       vel = velo_vs[i];
			}
			if (sdvn_vs[i] > 0.01) {
			   for (k=0;k<nsim_vs;k++) {
    			  dumv_vs = 0.0;
			      while (dumv_vs <= 0.0) {
					   qv_vs = random_vs();
					   xv_vs = invnor_vs(qv_vs);
					   dumv_vs = vel + xv_vs*sdvn_vs[i];
				  }
				  v2_vs[k] = dumv_vs;
			      vn_vs[k] = dumv_vs;
			   }
			}
			else {
			   for (k=0;k<nsim_vs;k++) {
				  v2_vs[k] = vel;
				  vn_vs[k] = vel;
			   }
			}

//           Calculate deltaz,depth and pseudo velocity			
			for (k=0; k<nsim_vs; k++) {
			   deltat_vs = t2_vs[k] - t1_vs[k];
			   deltaz_vs = calc_deltaz(ivm_vs[i],deltat_vs,vn_vs[k],vn_vs[k],
			                                     rka_vs[i],z1_vs[k]);
			   z2_vs[k] = z1_vs[k] + deltaz_vs;
			   v2_vs[k] = 2.0 * deltaz_vs/(t2_vs[k]-t1_vs[k]);
			   
    		}

//           Calculate mean and sd           
		    sum_vs  = 0.0;
		    sum2_vs = 0.0;
			double sumsqr = 0.0;
		    for (k=0;k<nsim_vs;k++) {
		       sum_vs  = sum_vs + z2_vs[k];
			   tsq_vs  = z2_vs[k];
		       sum2_vs = sum2_vs + tsq_vs*tsq_vs/(double) nsim_vs;                            
		    }
		    rav_vs[i+2] = sum_vs / (double) nsim_vs;
			for (k=0;k<nsim_vs;k++) {
				sumsqr = sumsqr + ((z2_vs[k] - (sum_vs/(double) nsim_vs))*(z2_vs[k] - (sum_vs/(double) nsim_vs)));
			}
			tsq_vs      = sumsqr / (double) nsim_vs;
            ssq_vs      = sum2_vs - tsq_vs*tsq_vs;
            if(ssq_vs<0.0) {
                rsd_vs[i+2]=0.0;
            }
            else {
			  rsd_vs[i+2] = sqrt(tsq_vs);
			}
         
		 } /* end of if(i>0&& deltat_vs < .....   else  */
	  }/* End of loop per Layer */

//     Write results to file
	  double xround,yround;
	  xround = (round(rav_vs[0]*10.0))/10.0;
	  yround = (round(rav_vs[1]*10.0))/10.0;
	  if (nfirst_vs != 0) {
         fprintf(fd_avz_vs,"\n%10.1f%10.1f",xround,yround);
      }
      else {
         fprintf(fd_avz_vs,"%10.1f%10.1f",xround,yround);
      }

	  for (k=0;k<nlay_vs;k++) {
		 sprintf(string3,"%8.1f",rav_vs[k+2]);
		 string3[strlen(string3)] = '\0';
		 if (strcmp(string3,"    NaN") == 0) {
			rav_vs[k+2]=9999.9;
		 }
		 fprintf(fd_avz_vs,"%8.1f",rav_vs[k+2]);
	  }
      xround = (round(rsd_vs[0]*10.0))/10.0;
	  yround = (round(rsd_vs[1]*10.0))/10.0;
	  if (nfirst_vs != 0) {
	     fprintf(fd_sdz_vs,"\n%10.1f%10.1f",xround,yround);
      }
      else {
	     fprintf(fd_sdz_vs,"%10.1f%10.1f",xround,yround);
      }
      nfirst_vs  = 1;

	  for (k=0;k<nlay_vs;k++) {
		 sprintf(string3,"%8.1f",rsd_vs[k+2]);
		 string3[strlen(string3)] = '\0';
		 if (strcmp(string3,"    NaN") == 0) {rsd_vs[k+2]=9999.9;}
		 fprintf(fd_sdz_vs,"%8.1f",rsd_vs[k+2]);
	  }
	} /* end of if (!skip_point) part)  */
  } /* end of loop per xy point */

//  Free memory, Close Files
  free(z1_vs);
  free(z2_vs);
  free(vn_vs);
  free(v1_vs);
  free(v2_vs);
  free(t1_vs);
  free(t2_vs);

  fclose(fd_grids);
                       
  fprintf(fd_avz_vs,"\n");
  fprintf(fd_sdz_vs,"\n");
  fclose(fd_avz_vs);
  fclose(fd_sdz_vs);

  fprintf(fd_rep_vs,"\nThe number of Points processed is : %4d\n",np_vs);
  fprintf(fd_rep_vs,"\nNORMAL EXIT SIMULATION\n");
  fprintf(fd_rep_vs, "Start Simulation %s %s\n", startdateStr, starttimeStr);
  char enddateStr [9];
  char endtimeStr [9];
  _strdate( enddateStr);
  _strtime( endtimeStr );
  fprintf(fd_rep_vs, "End Simulation   %s %s\n", enddateStr, endtimeStr);
		    
  fclose(fd_rep_vs);
  fclose(fp_log_vs);
             
   exit(0);
}

/********************************************************************
 *   Function : invnor_vs
 *
 *   Function to compute a function value (y * stdev)
 *   given a probability for a gaussian distribution
 *
 *   Derived from www.johndcook.com/normal_cdf_inverse.html 
 *   (based on Abrhamovic)
 *
 ********************************************************************/

double invnor_vs(double q)
{
   double tt,x,y;
   double c0 = 2.515517;
   double c1 = 0.802853;
   double c2 = 0.010328;
   double d1 = 1.432788;
   double d2 = 0.189269;
   double d3 = 0.001308;


   if (q <= 0.5) {
        tt = sqrt(-2.0*log(q));
   }
   else  {
        tt = sqrt(-2.0*log(1.0-q));
   }
   y = tt - ((c0 + ((c1 + c2*tt)*tt))/
	    (1 + (d1 + (d2 + d3*tt)*tt)*tt));

   if (q <= 0.5) {
      return(y);
   }
   else {
      return(-1.0*y);
   }
}

/********************************************************************
 *  Function random_vs
 *
 *  Function to generate a random number between 0.0 and 1.0
 *  
 ********************************************************************/
 
double random_vs()
{
    double x;
    int    i;
    x = 0.0;
    while (x < 0.00000001 || x > 0.999999999) { 
       x = (double) rand() / (double) RAND_MAX;
    }
    return x;
}

/********************************************************************
 *  Function calc_deltaz
 *
 *  Function to calculate delta z given delta t and velocity model
 *  parameters
 *  Input velmodel : velocity model of layer (1=Vint,2=Vreg,3=V0 and k)
 *                                            4=Vreg and k
 *        deltat   : delta two travel time time(s).
 *        velo     : fixed layer velocity(m/s)(model 1 and 3)
 *        vreg     : regeonally varying layer velocity (model 2 and 4)
 *        rka      : K (model 3 and 4)
 *        z1       : average depth top layer (m)
 *
 *  Output deltaz (m).
 *  Error condition : deltaz=-1.0 (case K =<0.0000001);
 *
 *  
 ********************************************************************/
double calc_deltaz(int velmodel, double deltat, double velo, double regv,
				                 double rka, double z1)
{
	double deltaz;
	double f,g;
    switch(velmodel)
    {
       case CONS_V: 
	       deltaz = deltat * velo /2;
	       break;

	   case REG_V: 
           deltaz = deltat * regv /2;
		   break;

	   case CONS_V_K:
		   if (rka>0.0000001) {
		       deltaz = ((z1*rka + velo)/rka)*
				  (exp( (double) (rka*deltat/2))-1);
		   }
		   else {
			   deltaz = -1.0;
		   }
           break;

	   case REG_V_CONS_K:
		   if (rka>0.0000001) {
			   deltaz = ((z1*rka + regv)/rka)*
				  (exp( (double) (rka*deltat/2))-1);
		   }
		   else {
			   deltaz = -1.0;
		   }
		   break;
           
	   case LIN_REGR:
		   f = velo;
		   g = rka;
		   deltaz = (f * deltat/2) + g;
           break;

	   default:
          break;
    }
    return deltaz;
}
double round(double d) 
{ 
  return floor(d + 0.5); 
} 